#include "classes.h"

#define EM_E0 GSL_CONST_MKSA_VACUUM_PERMITTIVITY //8.854187817e-12
#define EM_M0 GSL_CONST_MKSA_VACUUM_PERMEABILITY //1.25663706144e-6

double Structure::cosa(double fis)
{
    double val;
    val=sin(fi-fis)*sin(teta);
    return val;
}

double Structure::cosb(double fis)
{
    double val;
    val=sin(fi-fis)*cos(teta);
    return val;
}

double Structure::cosg(double fis)
{
    double val;
    val=cos(fi-fis)*sin(teta);
    return val;
}

double Structure::cosd(double fis)
{
    double val;
    val=cos(fi-fis)*cos(teta);
    return val;
}

__complex__ double Structure::ji1(double ro)
{
    __complex__ double val;
    val=1I*gsl_sf_bessel_I1(__imag__(alfa2)*ro);
    return val;
}

__complex__ double Structure::yi1(double ro)
{
    __complex__ double val;
    val=-1*gsl_sf_bessel_I1(__imag__(alfa2)*ro)+2I/M_PI*gsl_sf_bessel_K1(__imag__(alfa2)*ro);
    return val;
}

double Structure::js1(double ro)
{
    double val;
    val=gsl_sf_bessel_J0(alfa1*ro)-gsl_sf_bessel_J1(alfa1*ro)/alfa1/ro;
    return val;
}

double Structure::ys1(double ro)
{
    double val;
    val=gsl_sf_bessel_Y0(alfa1*ro)-gsl_sf_bessel_Y1(alfa1*ro)/alfa1/ro;
    return val;
}

double Structure::js2(double ro)
{
    double val;
    val=gsl_sf_bessel_I0(__imag__(alfa2)*ro)-gsl_sf_bessel_I1(__imag__(alfa2)*ro)/__imag__(alfa2)/ro;
    return val;
}

__complex__ double Structure::ys2(double ro)
{
    __complex__ double val,y0;
    y0=1I*gsl_sf_bessel_I0(__imag__(alfa2)*ro)-2/M_PI*gsl_sf_bessel_K0(__imag__(alfa2)*ro);
    val=y0-yi1(ro)/alfa2/ro;
    return val;
}

__complex__ double Structure::cexp(double rr)
{
    __complex__ double val=1I;
    val=cos(k*rr)-val*sin(k*rr);
    val/=rr;
    return val;
}

__complex__ double Structure::cexp0(double ro, double fis)
{
    __complex__ double val=1I;
    val=cos(k*ro*cosg(fis))+val*sin(k*ro*cosg(fis));
    return val;
}

__complex__ double Structure::G(int n)
{
    __complex__ double val;
    float ee;
    if(n==1) ee=e1; else ee=e2;
    val=ee*js1(a)/gsl_sf_bessel_J1(alfa1*a)/alfa1-Rs1(a)/R1(a)/alfa2;
    val*=w*EM_E0*a/beta/(1/pow(alfa1,2)+1/pow(__imag__(alfa2),2));	//Только для чистомнимых alfa2
    return val;
}

double Structure::R1(double ro)
{
    double val;
    //val=ji1(ro)*yi1(b)-ji1(b)*yi1(ro);
    I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
    K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
    val=2/M_PI*(I1b*K1r-I1r*K1b);
    return val;
}

__complex__ double Structure::Rs1(double ro)
{
    __complex__ double val;
    //val=js2(ro)*yi1(b)-ji1(b)*ys2(ro);
    I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
    K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
    I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
    K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
    ro=__imag__(alfa2)*ro;
    val=2I/M_PI*(I0r*K1b-I1r*K1b/ro+I1b*K0r+I1b*K1r/ro);
    return val;
}

__complex__ double Structure::S1(double ro)
{
    __complex__ double val;
    //val=ji1(ro)*ys2(b)-js2(b)*yi1(ro);
    I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
    K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
    I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
    K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
    ro=__imag__(alfa2)*b;
    val=2I/M_PI*(I1r*K1b/ro-I1r*K0b+I1b*K1r/ro-I0b*K1r);
    return val;
}

double Structure::Ss1(double ro)
{
    double val;
    double rob;
    rob=__imag__(alfa2)*b;
    //val=js2(ro)*ys2(b)-js2(b)*ys2(ro);
    I1r=gsl_sf_bessel_I1(__imag__(alfa2)*ro);
    K1r=gsl_sf_bessel_K1(__imag__(alfa2)*ro);
    I0r=gsl_sf_bessel_I0(__imag__(alfa2)*ro);
    K0r=gsl_sf_bessel_K0(__imag__(alfa2)*ro);
    ro=__imag__(alfa2)*ro;
    val=2/M_PI*((I0b-I1b/rob)*(K0r+K1r/ro)-(I0r-I1r/ro)*(K0b+K1b/rob));
    return val;
}

__complex__ double Structure::Ero1(double ro,double fis)
{
    __complex__ double val;
    val=1I*(-beta*alfa1*js1(ro)+w*EM_M0*G1*gsl_sf_bessel_J1(alfa1*ro)/ro)*sin(fis);
    return val;
}

__complex__ double Structure::Efi1(double ro, double fis)
{
    __complex__ double val;
    val=1I*(-beta*gsl_sf_bessel_J1(alfa1*ro)/ro+w*EM_M0*G1*alfa1*js1(ro))*cos(fis);
    return val;
}

__complex__ double Structure::Hro1(double ro, double fis)
{
    __complex__ double val;
    val=1I*(-beta*G1*alfa1*js1(ro)+w*EM_E0*e1*gsl_sf_bessel_J1(alfa1*ro)/ro)*cos(fis);
    return val;
}

__complex__ double Structure::Hfi1(double ro, double fis)
{
    __complex__ double val;
    val=1I*(beta*G1*gsl_sf_bessel_J1(alfa1*ro)/ro-w*EM_E0*e1*alfa1*js1(ro))*sin(fis);
    return val;
}

__complex__ double Structure::Ero2(double ro,double fis)
{
    __complex__ double val;
    val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta/R1(a)*Rs1(ro)+w*EM_M0*G2/ro/alfa2/S1(a)*S1(ro))*sin(fis);
    return val;
}

__complex__ double Structure::Efi2(double ro, double fis)
{
    __complex__ double val;
    val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta/R1(a)/alfa2*R1(ro)/ro+w*EM_M0*G2/S1(a)*Ss1(ro))*cos(fis);
    return val;
}

__complex__ double Structure::Hro2(double ro, double fis)
{
    __complex__ double val;
    val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(-beta*G2/S1(a)*Ss1(ro)+w*EM_E0*e2/alfa2/R1(a)*R1(ro)/ro)*cos(fis);
    return val;
}

__complex__ double Structure::Hfi2(double ro, double fis)
{
    __complex__ double val;
    val=1I*pow(alfa1,2)/alfa2*gsl_sf_bessel_J1(alfa1*a)*(beta*G2/alfa2/S1(a)*S1(ro)/ro-w*EM_E0*e2/R1(a)*Rs1(ro))*sin(fis);
    return val;
}

double Structure::alfa()
{
    int j=0,q;
    double alf=1000,xa,xb,det,res=0;
    __complex__ double d1,d2;
    det=10;
    while (j<1)
    {
        xa=alf;	xb=xa+det;
        //first point
        alfa1=xa;
        beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
        alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
        alfa2=1I*alfa2;
        I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
        I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
        K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
        K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
        d1=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
        d1=d1*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
        d1=d1-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
        //second point
        alfa1=xb;
        beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
        alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
        alfa2=1I*alfa2;
        I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
        I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
        K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
        K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
        d2=alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a);
        d2=d2*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
        d2=d2-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
        if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))
        {
            for (q=0;q<16;q++)
            {
                alf=(xa+xb)/2;
                // first point
                alfa1=xa;
                beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
                alfa2=1I*alfa2;
                I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
                I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
                K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
                K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
                d1=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
                d1=d1*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
                d1=d1-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
                //second point
                alfa1=alf;
                beta=sqrt(pow(k,2)*e1*m1-pow(alfa1,2));
                alfa2=sqrt(pow(beta,2)-pow(k,2)*e2*m2)+0I;
                alfa2=1I*alfa2;
                I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
                I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
                K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
                K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
                d2=(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a));
                d2=d2*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a));
                d2=d2-pow(beta/k,2)*pow(1+pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1;
                if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))	xb=alf;
                else	xa=alf;
            }
            res=alf;
            j++;
        }//j++;
        alf=alf+det;
    }
    //return __real__((alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-m2/m1*pow(alfa1*a,2)/alfa2/a*Ss1(a)/S1(a))*(alfa1*a*js1(a)/gsl_sf_bessel_J1(alfa1*a)-e2/e1*pow(alfa1*a,2)/alfa2/a*Rs1(a)/R1(a))-(beta/k)*pow(1-pow(alfa1,2)/pow(__imag__(alfa2),2),2)/e1/m1);
    return	res; //__real__(d1/d2);
    //return __imag__(S1(a));
}
